
#include <iostream>

#include "polynomial_base.hh"
#include <boost/multi_array.hpp>

int main ()
{
  // Define structure of tensor
  const int dim = 3; // 3d matrix
  typedef int Index;
  typedef boost::multi_array<double, dim> Tensor;
  Tensor tensor;
  
  // Define the extents of each separate dimension
  const int NT = 3;
  const int NSIM = 4;
  const int NDIM = 3;
  tensor.resize(boost::extents[NT][NSIM][NDIM]);

  // Accessing elements using []
  for (Index i = 0, offset; i != NT; ++i)
  {
    for (Index j = 0; j != NSIM; ++j)
    {
      for (Index k = 0; k != NDIM; ++k, ++offset)
      {
        tensor[i][j][k] = offset;
        std::cout << "Tensor[" << i << "][" << j << "][" << k << "] = " << tensor[i][j][k] << "\n" ;
      }
    }
  }

  return 0 ;
}
